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^ ■ ABSTRACT 

I I There is currently no accepted theoretical framework for the formation of the most 

O . massive stars, and the manner in which protostars continue to accrete and grow in mass 

beyond ^lOM© is still a controversial topic. In this study we use several prescriptions 
^ I of stellar accretion and a description of the Galactic gas distribution to simulate the 

luminosities and spatial distribution of massive protostellar population of the Galaxy. 
We then compare the observables of each simulation to the results of the Red MSX 
Source (RMS) survey, a recently compiled database of massive young stellar objects. 
• We find that the observations are best matched by accretion rates which increase 

\l I as the protostar grows in mass, such as those predicted by the turbulent core and 

. competitive accretion (i.e. Bondi-Hoyle) models. These 'accelerating accretion' models 

' provide very good qualitative and quantitative fits to the data, though we are unable 

. I to distinguish between these two models on our simulations alone. We rule out models 

■ with accretion rates which are constant with time, and those which are initially very 

' high and which fall away with time, as these produce results which are quantitatively 

and/or qualitatively incompatible with the observations. To simultaneously match the 
low- and high-luminosity YSO distribution we require the inclusion of a 'swollen-star' 
pre-main-sequence phase, the length of which is well-described by the Kelvin-Helmholz 
timcscale. Our results suggest that the lifetime of the YSO phase is ^ lO'^yrs, whereas 
^ the compact Hii-region phase lasts between 2 — 4 x lO^yrs depending on the final 

I mass of the star. Finally, the absolute numbers of YSOs are best matched by a globally 

averaged star-formation rate for the Galaxy of 1.5-2AfQyr~^. 

Key words: stars: massive, stars: pre-main-sequence, stars: protostars, stars: forma- 
tion 



1 INTRODUCTION 

The formation mechanism for stars with masses >10Mn 
is still unclear. The steepness of the Initial Mass Func- 
tion (IMF) means that such stars are rare, while the short 
Kelvin-Helmholz contraction times of massive stars mean 
that they arrive on the main-sequence (MS) whilst still ac- 
creting and heavily embedded within their natal molecular 
clouds. Meanwhile, from a theoretical point-of-view there is 
considerable doubt that the 'classical' theory of star for- 
mation is applicable massive stars. Firstly, in the model 
of an isothermal collapsing sphere, the ac cretion rat e de- 
pends only on the sound-speed of the gas l|Shulll977l ). and 



implies unrealistically large formation tim escales for stars 
with masses >1OM0 (|Stahler et al.l I2OO0I ). Secondly, once 
the massive proto-star arrives on the MS it will exert con- 
siderable outward radiation pressure on the infalling ma- 
terial. In the spherically-symmetric case, this can impede 
and eventually halt accretion, preventing the formation of a 
high- mass star (|Kahnlll974l : IWolfire fc Cassinellilll987l ). For 
further discussions regarding both the theoretical and ob- 
servational challenges faced in the stu dy of massive star for- 
mation, see IZinnecker fc Yorkd (|2007l ). 

There are several proposed amendments to the classi- 
cal model of star formation to account for massiv e stars, a 
few of which we now discuss. iMcKee fc Taiil (|2003l . hereafter 
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MT03) suggested a mechanism whereby the larger molecular 
pre-stellar cores are supported against collapse by larger tur- 
bulent motions (the turbulent core model) . Once these cores 
become unstable to collapse they have accretion rates which 
accelerate with time. In the competitive accretion model 
l|BonneU fc Batel2006l ) , a large reservoir of gas is funnelled to 
the centre of a cluster of protostars, allowing the cores at the 
cluster's centre to achieve ve ry large accretion ra t es thr ough 
the Bondi-Hoyle mechanism. ISchmeia fc KlessenI (|2004l ') pro- 
posed that the accretion rates of proto-stars within a gas 
cloud are governed by the complex interplay of self-gravity 
and turbulence, and that upon collapse they have initially 
very high accretion rates owing to the large amounts of ma- 
terial that has accumulated in their immediate viscinity (the 
gravo-turbulent model). Indeed, several authors have pro- 
posed models for star formation in which accretion rates 
behave in this way with time JPostc r & Cho valieil 1 19931: 
Whitworth fc Ward-Thompson|[200ll: iMotovama fc Yoshidal 



20031 ). Meanwhile, 3-D hydrodynamical simulations of mas- 
sive forming stars through accretion from a circumstellar 
disk displayed accretion rates which were approximately 
constant throughout the simulation, though they may pass 
through distinct accretion phases and be highly variable o n 
short timescales (|Krumholz et al.ll2009l : iKuiper et al.ll20lil ). 

The mechanisms described above predict accretion rates 
which evolve differently throughout the formation timescale 
of the star. The turbulent core and competitive accretion 
(Bondi-Hoyle) models, though they differ in form, both pre- 
dict accretion rates M that should accelerate with time as 
the central star grows in mass. The gravo-turbulent model 
however predicts that there should be an initial burst of ac- 
cretion, which gradually tails off as the star approaches its fi- 
nal mass. These can be contrasted to the 'standard model' of 
star-forma tion, as well as th ose found in the numerical sim - 
ulations of lKrumholz et al.1 (2009) and lKuiper et al.l (|201ll ). 
when the accretion rate is largely constant. As the accretion 
rate and its evolution will govern the proto-star's luminos- 
ity as a function of time, as well as the star's formation 
timescale, with a large complete sample of massive form- 
ing stars it should be possible to distinguish between these 
scenarios from analysis of their luminosity distribution. 

1.1 The RMS Survey 

Such a sample has recently come to fruition. The Red 
MSX Source (RMS) survey l|Hoare et al.l l2005lfl has se- 
lected point-sources from the MSX catalogue fe gan et al.l 
l2003f) with colours appropriate for young massive stars 
i Lumsden et al.ll2002l ). and weeded out contaminant objects 
such as planetary nebulae and evolved stars through follow- 
up near- infrared spectroscopy and con tinuum radio obser- 
vations (|Urauhart et al.ll2007al . |2009b). Sources have kine- 
matic distance estimates from molecu lar line observations 
l|Urauhart et al.1 l2007bl . |2008| , l2009al ). while the sources' 
bolometric luminosities have been determined from fits t o 
their spectral energy distributions (jMottram et al.ll2011bl ). 
High-resolution mid-infrar ed imaging has been u sed to inves- 
tigate source multiplicity (|Mottram et al.ll2007l '). The RMS 
databcise now contains over 1,000 confirmed massive young 
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stars. These objects are separated into massive Young Stel- 
lar Objects (YSOs) and compact Hll-regions based on one 
of two criteria. Objects which are spatially extended in the 
high resolution mid-infrared images (spatial resolution ~l") 
due to Lyman-a heating of the surrounding dust, or which 
have detectable radio emission (>lmjy/beam) are consid- 
ered to be ionizing their surroundings, and hence are clas- 
sified as Hll-regions. These objects are considered to have 
reached the main-sequence, though they may still be accret- 
ing. Objects displaying no evidence for an ionized nebula are 
classified as YSOs. 

In this paper we use the results of the RMS survey to 
test the predictions of different accretion laws. We do this 
by constructing realistic populations of massive protostars 
and distributing them throughout the Galaxy according to 
a model of the Galaxy's gas distribution. We compute the 
observed properties of the objects in the simulation, apply 
the RMS selection and classification criteria and compare 
them to the observed results. 

This pape r is similar in concept to that of 
[Froebrich ct al.' (2006'), who used the numerical results 
of Schmoia & Klcsscn (2004) to make predictions about 
the earliest stages of star f ormation. It is also similar to 
iRobitaille fc Whitnevl (|201G| ). who calculated the evolution 
of the spectral energy distribution (SED) of intermedi- 
ate mass YSOs M^'^S-IOMq using a single accretion law. 
From this, they calculated the predicted fiuxes in the four 
Spitzer/IRAC channels and compared the total number of 
objects to those found in the GLIMPSE survey. They then 
used this information to predict the global star-formation 
rate of the Galaxy. Our approach is different, in that we 
predict the 21/^m luminosity function of massive YSOs for 
a variety of different accretion scenarios, with the ultimate 
goal of understanding how massive stars form, although we 
too are able to put some constraints on the Galactic global 
star-formation rate. 

We begin in Sect. [2] with a description of the model 
construction. The results of the model, the effect of the free 
parameters, and of the different accretion laws, are presented 
in Sect. O We discuss the results and the implications for 
models of massive star formation in Sect. |4] We provide a 
summary of our results in Sect. [5] 



2 MODEL CONSTRUCTION 

The methodology for the construction of the model is as fol- 
lows: assuming a global star formation rate of the Galaxy, we 
generate lO^yrs of stars - up until the point where the oldest 
massive stars have comfortably reached the main-sequence 
and have broken free of their natal environments. Using a 
description of the gas distribution in the Galaxy, these stars 
are each inserted into the model Galaxy at a random loca- 
tion dictated by the relation between local gas density and 
local star-formation rate. That is, locations of high den- 
sity have higher star formation rates, and hence are more 
likely to host newly-formed stars. Quantitatively, we scale 
the star-formation rate with gas density SFRoc , in ac- 
cordance with_th«_Sclinu law for star-forming 
galaxies l|Kennicuttl 1 19981 ) . Typically, the value of « 1.4 
is quoted, though some authors have suggeste d that A^ is 
closer to unity, and may even be as low as 0.8 (|Blanc et al.l 
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I2009I ). Strictly this law applies to the surface density Edens 
of galaxies, but since the scale height of the Galaxy's gas 
distribution is small compared to the radial size, p oc Edens- 
Each star in the simulation is given an age (i.e. time 
since accretion first began), which is randomly selected from 
a uniform distribution between - lO^yrs. Using a given 
model of stellar accretion, the age, and final mass Man of 
each star, we compute the current mgiss Mcur - The pre-main- 
sequence birthlines of iHosokawa fc Omukail (|2009l ) are then 
used to compute the luminosities of each star, which we then 
use to calculate the observed flux at 21^m. For stars which 
we determine have arrived on the main-sequence, making 
standard assumptions we calculate the properties of the sur- 
rounding Hll-region, and ultimately the observed radio flux. 
We then pass the observables through the RMS selection 
criteria, and those sources which would be detected in the 
survey are then classified according to their observed prop- 
erties. Below, each step of the process is described in more 
detail. 

2.1 Galactic gas distribution 

For a prescription of the distribution of baryonic material 
in the Galaxy, we use as a basi s the Galact ic distribution of 
free electrons as determined bv lTavlor fc Co rdcs ( 199 3, and 
refs therein), and later updated by ICord^ fc Lazior(|2002l ) 
in a model known as NE2001. These studies use pulsars 
with known distances to infer the line-of-sight density of free 
electrons. We assume that in the Milky Way disk, the total 
ISM density (neutral gas, dust etc.) scales directly with the 
density of free electrons. 

The model consists of three components: thin disk, thick 
disk, and four spiral arms. For the thin and thick disk com- 
ponents we use the definitions and physical parameters as 
defined in NE2001. For the spiral arms, we have made slight 
adjustments: we used the coordinates of the fiducial points 
for each of the four arms as quoted by TC93 and fitted each 
arm with a logarithmic spiral, 

= Aspj log(rj/r„i„,j) -I- eo,j (1) 

where 9j and rj define the Galactocentric distance and az- 
imuthal anglqj of each spiral 9oj are fitted 

constants; and rminj is the minimum Galactocentric radius 
of the defining fiducial points in TC93. We specified that 
two spiral arms should originate from the ends of the cen- 
tral Bar. To ensure this we adjusted the fiducial points such 
that Vmin,] was always 4kpc (the length of the Bar as de- 
termined by .Benjamin et al. 200 5) and 6*0, j was 45° for the 
Scutum-Crux arm and 225° for the opposing (Perseus) arm. 
The other two arms (Sagittarius and Norma) were traced 
inwards according to the logarithmic fit, and the inner re- 
gions of the arms were attenuated to blend with the inner 
disk component using an attenuation function 5", similar to 
that used in the NE2001 model; 

5" = Sech^(r - dinner) for (r - Ainnor) < 
S= 1 for (r - Ainnor) > 

We used a somewhat arbitrary value of 5kpc for ylinncr in 
^ 8j =0 is defined as the Sun - Galactic Centre axis. 
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Figure 1. Gas surface densities of the model Galaxy in the XY 
(top) and XZ (bottom) planes. The panels show the linear grey- 
scale of the normalised gas density, with the relevant scale shown 
above each panel. The thick-disk, thin-disk and spiral-arm com- 
ponents, as described in TC93, are included. The location of the 
sun is indicated on the top panel. 

order to blend them with the inner thin disk. For the outer 
disk, the spiral arms were extrapolated radially outwards 
to 21kpc, at which point their contribution to the local gas 
density is negligible. The normalised surface density of the 
Galaxy gas distribution model in both the XY and XZ 
planes can be seen in Fig. [T] 

The description of the Galactic gas density we use has 
not been precisely tuned to the current state-of-the-art pic- 
ture of our Galaxy's morphology. For example, we do not 
include the 'kinks' in the arms in the Solar neighbourhood, 
nor do we include other smaller features such as the 'local' 
arm. Our goal here is not to accurately predict the spatial 
distribution of the YSO population, nor is it to trace the 
spiral structure of our Galaxy. The purpose of this work 
is to model the observed luminosity distributions of proto- 
stars, and so we only need a description of the Galaxy's 
morphology which is realistic in terms of size-scales (e.g. 
the scale-height, z) and discrete structures (e.g. the num- 
ber and size-scales of spiral arms). Tweaks that we could 
make to the model gas distribution to better match the 
spatial distribution will be discussed later, while a better 
match to the 3-dimensional distribution of young massive 
stars aw aits the impleme ntation of more accurate distances 
(see e.g. iReid et alflioogl ). 



4 B. Davits et al. 




log T eff log T eff 

Figure 2. Evolution tracks of massive pre-main-sequence stars for the original accretion rates of M = lO~'*M0yr^^ (upper panel) and 
lO~'^M0yr~^ (lower panel) in the H-R diagram. The dotted line represents the birth line at these accretion rates in each panel. In the 
upper (lower) panel, the solid lines present the evolution tracks of 12, 10, 8, and 6Mq (25, 20, 15, and IOMq) stars after mass accretion 
ceases. The symbols on each track indicate the elapsed time since the ter mination of mas s accretion, 10*, 3 X lO'*, 10^, and 3 X 10^ years. 
The dashed line presents locii of zero-age main-sequence stars taken from lSchaller et alj l l 19921 ). Stellar radius is constant at -R* = 1 Rq, 
10 -Rq, and 100 Rq along the dotted lines. 



2.2 Generation of stellar populations 

Assuming a global star-formation rate for the Galaxy of 
SFRcai, we generate IMyr of star- formation, under the as- 
sumption that any star older than IMyr would be sufficiently 
free of its natal material to be precluded by the RMS selec- 
tion criteria as being too evolved. Estimates of SFRgs_i range 
from ~2-lOM0yr~^, though m ost measure ments appear to 
be in the range 3-4M0yr~^ (see lPiehl et al.ii,2006, . and refer- 
ences therein). Rec ently, in a study which was si milar to the 
one presented here, iRobitaille fc WhitnevI ()2010l ) attempted 
to derive SFRgs_i from a sample of YSOs selected from 
GLIMPSE, finding that values of S'™Gai=O.7-1.5M0yr-^ 
were preferred. This is somewhat lo wer than the other es- 
timate s described above. However, IRobitaille fc WhitnevI 
(|201[1 ) used only one description of accre tion h istory in 
their study, that of |Bernascoiii^ Macde^ 1 19961 ). though 
it is unclear whether they used the constant accretion 
rates or those that i ncreased with protostellar mass (see 
iRobitaiiie" et al.ll2006l'). In eithe r case, these accretion rates 
of Bernasconi fc Maedeij l| 19961 ) are an order of magnitude 
lower than those presented in more contemporary studies of 
e.g. MT03. This produces slightly longer pre-MS lifetimes 
and therefore results in larger numbers of YSOs. This may 
explain in part why IRobitaille fc WhitnevI (|201Cll ) required 
values of SFRq^i which were lower than in previous stud- 
ies. For now we adopt 5'-Fi?Gai=3M0yr~^. This parameter 
affects only the overall normalization of objects, and the ef- 
fect on our conclusions of allowing this parameter to vary 
slightly will be explored later. 

For the total mass of stars in our model, we generate a 
distribution of stel lar mass e s acco rding to the Initial Mass 
Function (IMF) of lKroupal l|200ll ). which is essentially the 
Salpeter law but augmented at low masses. The masses we 



generate are the final stellar masses M^n, i.e. the mass that 
each star will have once it has finished accreting. We make 
the assumption that the global star formation rate has been 
constant over the last lO^yrs, and the age of each star in the 
simulation is selected from a uniform random distribution of 
ages from zero to lO^yrs. We did experiment with 'cluster- 
ing' the star formation into localized bursts, which is likely 
to be a more realistic representation of how star formation 
proceeds within the Galaxy, but this was found to have no 
impact on the results other than to increase stochastic ef- 
fects within a single simulation. 



2.3 Computation of observables 

In order to calculate the observed properties of the 
protostars in our simulations, we utilize the theoreti- 
cal evolution tracks for a ccreting protostars taken from 
iHosokawa fc Omukail (|2009l . hereafter HO09). The authors 
calculated the evolution of accreting protostars for various 
constant accretion rates. They show that the evolution with 
high accretion rates of > 10~* Mq yr~^ has some charac- 
teristic features. With 10~^ Mq yr~^, for example, the pro- 
tostar becomes significantly swollen and does not reach the 
main-sequence (MS) for M* < Mms ^ SOMq. The critical 
mass AfMS becomes higher with the higher accretion rate. If 
the mass accretion ceases before the arrival to the MS, the 
star experiences a contraction phase lasting approximately 
a Kelvin-Helmholtz timescale tKU, 



GMi 



(2) 



This predicts that massive pre-main-sequence (PMS) 
stars exceeding 8 Mq should exist. Figure [2] shows the evo- 
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lution tracks for such massive PMS stars in the HR dia- 
gram. In this paper, we assume that the massive PMS stars 
contract keeping the constant himinosity for simphfying the 
detailed evolution. 

2.3.1 Accretion laws 

To calculate the physical properties of the forming star, we 
must first know how much matter has been accreted since 
the star began to form. For this, we require a prescription 
of the star's accretion rate. In this study we investigate 
three distinct accretion-rate scenarios: constant accretion- 
rate; accretion-rates which increase as the mass of the star 
grows; and those which decrease. Below, we describe the pre- 
scriptions for each of these accretion rates, while in Fig. [3] 
we illustrate the behaviour of the accretion rates with both 
mass and time for each of these laws. 

a) Constant accretion Here we use the results of IIO09 
directly. For the age t and final mass Mfln of each star in 
our simulation, we calculate the current mass for a specified 
accretion rate M from, 

Afcur = Mt (3) 

We can also calculate the formation timescale of the star 
tform from, 

tform = Mfin/M (4) 

We use the IIO09 tracks to calculate the luminosity and 
temperature of each star from its current mass Mem. For 
stars which have reached their final mass but whose masses 
are below the Mms for that accretion rate, we calculate the 
Kelvin-Helmholz timescale txH. If a star is not older than 
(tform +tKH) thcu it is cousidcrcd to still be 'swollen', that 
is it has not yet contracted to the main-sequence, and so 
emits a negligible amount of ionizing fiux. Stars which have 
arrived on the main-sequence are considered to have Lyman 
fluxes similar to zero-age main-sequence stars of the same 
mass, but with extra luminosity due to accretion. 

b) Accelerating accretion I: MT03 In the MT03 pre- 
scription of accretion by a massive protostar in the turbulent 
core model, the accretion rate M and formation timescale 
iform of an accreting star are given by their Equs. (41) and 
(44), 

M^4.6xl0-Meyr-(^)°%r(^)°'^ (5) 

t... = 1.29xl0V(^)""E-- (6) 

where Ed is the surface density of a prestellar clump in units 
of gcm~^. From these equations we can find an expression 
for the current mass Mcur of an accreting star, 

Using these equations, and the randomly generated stellar 
ages, we can then calculate the current mass and formation 
timescale of each star. 



To calculate the luminosities of the accreting star, we 
must take into account the stellar luminosity and the accre- 
tion luminosity. To do this we take the star's current mass 
and current accretion rate, and look up the corresponding 
luminosity from the tracks of HO09. In doing so, we are 
making two approximations: firstly, we linearly interpolate 
the HO09 tracks in log-space at O.Oldex increments between 
10~^ and 10~^MQyr~^, since IIO09 only calculate tracks at 
10"^ 10"* and W'^Meyi-'^. Secondly, we are assumme 
that a star whose accretion rate has accelerated to the cur- 
rent value has a similar luminosity to a star whose accretion 
rate has always been constant at the current value. While 
one may expect there to be systematic differences between 
these two cases, we found that the Afcm-J^boi relationship of 
our approximation was quantitatively similar to that com- 
puted by MT03 (their Fig. 6) , since the accretion luminosity 
only contributes a small fraction to the total luminosity. 
For stars which have finished accretin g, luminosties are 
taken from the stellar structure models of lMevnet fc Maederl 
(|2000l ). 

Again, MT03 find that very massive stars join the MS 
at Mms ~20Mq, though the precise value of Mms depends 
on Eel. As with the case of constant accretion, stars which 
have finished accreting but have final masses below Mms 
are assumed to undergo a pre-MS contraction phase which 
approximates to tKH for that particular star. 

c) Accelerating accretion II: Bondi-Hoyle accretion 

Similar accelerating accretion is predicted under the Bondi- 
Hoyle (BH) accretion releva nt to the competitive accretion 
model PonneU fc Batj2006l '). The difference here is that the 
accretion rate depends only on t and Mcur, not on Mfin. The 
BH ac cretion rate is given by Equ. (6) of iBonnell fc Bate! 

where Vi and p are the velocity dispersion and gas density 
within the star-forming clump, and Mi is the initial mass of 
a star-forming clump that becomes unstable to gravitational 
collapse. Integrating and rea rranging this equation, and us- 
ing the values determined bv lBonnell fc Batd (|2006l ') for the 
parameters of Vi, Mi and p (0.5km s~^, 10~^^gcm'^, and 
0.5Mq respectively), the current mass of a star accreting in 
this way can be shown to be given by, 

Mcur = O.5M0 (t/loVr + 1)^ (9) 
The formation timescale is then simply, 

^--=^°'>'^ 

There has yet to be a thorough computation of BH birth- 
lines, so we must make assumptions for the object's luminos- 
ity, the duration of the pre-MS 'swollen' phase and the mass 
at which very massive stars arrive on the MS, Mms- The 
luminosity of a BH accreting object is again assumed to be 
well represented by the HO09 track with the same current 
accretion rate at the same Mcur. We again assume that the 
MS contraction time is equal to tKH, and we leave Mms as 
a free parameter in this model. 
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Figure 3. Evolution of the accretion rate for star with a final mass M^n of SOMq for the different prescriptions of accretion studied 
here. The panels show the accretion rate as a function of current mass Mcur (left) and as a function of time (right). The MT03 models 
plotted use St.i=lgcm^. 



d) Decelerating accretion For a representitive model of 
decelerating acc retion, we use the 'g r avo-t urbulent fragmen- 
tation' model of[S^meir&!^liiii3 HqqI, hereafter SK04). 
In this model, protostars at the centre of a cluster's poten- 
tial well are subject to a greater amount of infalling ma- 
terial, and the initial accretion rates are very high. As the 
reservoir of material is depleted, the accretion-rate gradu- 
ally decreases, eventually falling to zero. Before preceeding 
with a description of how we incorporate this accretion sce- 
nario into our simulations, we first note that in SK04 no 
stars more massive that lOM© are formed, and that these 
authors' goals were not necessarily to simulate massive star 
formation. Nevertheless, we include it in our study as a tool 
with which to examine how decelerating accretion affects the 
relative numbers of massive protostars. 

For a prescription of such accretion, we use the recipes 
of SK04, which are taken from empirical fits to their numer- 
ical simulations. Here the accretion rate is given by, 

log(M) = (logMo)-te"'^" (11) 

T 

where Mq and r are constants. The constant Mq specifies 
the maximum accretion rate, which is reached after a brief 
period of rapid acceleration in M . Once the accretion rate 
reaches Mo, it falls away gradually on a timescale dictated 
by r. 

As a representative case of decelerating accretion, we 
use the simulation 'G2' from SK04, which uses a gaussian re- 
alization of the initial turbulence within the molecular cloud. 
Here, the value of logMo was found to be a linear function 
of logMfln, up to A/finS^lOM0. Specifically, from fits to their 
data, 

log(Mo) = 0.75 log A/fi„ - 4.7 (12) 

Though SK04 formed no stars with masses greater than 
IOMq, for our simulations we make the simple assumption 
that this relation can be extrapolated up to the most mas- 
sive stars in our simulation, Affln^lSOM©. The parameter r 



also appears to be a linear function of logA/fln, though there 
is considerable dispersion in this relation. From our linear 
fits to the results of SK04, we find, 

r/loVrs = Ti log A/fln + ro (13) 
with, 

ro = (4.6 ± 0.3) X lO^yrs 
n = (4.3 ± 2.0) X lO^yrs 

Therefore, for a star of given Affln and age t, Eqs. (|ll|l . (|12|l 
and (|13p can be used to calculate that star's accretion rate 
as a function of time since t=0. Integration of this function 
then yields the star's current mass Afcur. 

Using these parameters however, the formation times of 
massive stars become unrealisticly long, and no stars with 
Affln>15AfQ can be formed within lO'^yrs. Massive stars can 
be formed within a reasonable time, but only if the constants 
in Eq. (|13p are pushed to their extreme limits, such that the 
largest possible values of r result. This will be discussed 
further in Sect. |4] 

As with the other accretion laws studied, luminosities 
are determined from the pre-MS tracks of HO09, Mms is left 
as a free parameter, and the length of the pre-MS 'swollen' 
phase is assumed to be given by the KH timescale. 



2.3.2 21jj,m lumtnosity 

Having calculated Lboi from the accretion law, we can now 
determine the observed fiux through the MSX 21/im filter, 
on e criterion used for sou rce selection in the RMS survey. 
In lMottram et al.l (|2011bl ) it was found that the 21^m fiux 
F2\ of objects in the sample scaled as a constant fraction of 
the total bolometric flux Ftot according to, 

log = 1.43 ±0.27 (14) 

V 1'21 J 
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No significant correlation was found between Ftot/F2i{= 
C21) and any observable parameter in the RMS survey, such 
as IR colour or F2i. 

To convert our Lboi determined in in Sect. l2.3Tn to F21, 
we use the following relation, 

Lhoi A21 



F2 



■10 



(-A21/2.5) 



(15) 



where C21 is randomly drawn from a normal distribution 
of numbers with a mean of 1.43 and standard deviation of 
0.27 (following Ea. ll4|l : A21 is the interstellar extinction (see 
below), and A21 is the effective bandwidth of the MSX 21fj,m 
mter. 

As yet we have not dealt exlicitly with the issue of line- 
of-sight extinction. There are likely two sources of extinction 
to each source - interstellar, and circumstellar. Circumstel- 
lar extinction will vary largely from source to source, since 
it will depend on random properties such as cavity open- 
ing angle and inclination to line-of-sight. We have assumed 
that the factor C21 incorporates the attenuating effect of 
circumstellar extinction. 

In terms of interstellar extinction at 21/im towards each 
object, ^21, one may expect a systematic trend of this factor 
with object distance. We have made a detailed study of this 
factor and how it may impact our results, using Galactic star 
clusters as test points (see Appendix A). Though we do ac- 
count for this effect by incorporating a factor of lo'^'^^i/z.s) 
into Eq. 1151 we expect at most that this plays a very minor 
role in the results of our simulation, since A21 < 1.5 for even 
the most distant objects. 

The current model does not account for a mid-IR 
faint or 'dark' phase early on in the evolution of massive 
stars. Examples of lumino us young objects that are mid-IR 
faint/dark clearly do exist. lEllingsenI (|2006t ) found that 30% 
of methanol maser sources, that are invariably associated 
with massive YSOs, do not have counterparts in the 4-8 mi- 
cron GLIMPSE survey. Also, about 20% of the active cores 
in infrared da rk clouds (IRDCs) have a luminosity greater 
than 50001/0 (|Rathborne et al.ll20ldl l. Though we are not 
yet explicitly accounting for an IR 'dark' phase here, during 
the very early stages of the growth of massive stars their 
luminosity will be below 50001/0 and therefore not in the 
RMS survey; an aspect that is already accounted for in the 
modelling presented here. Since very few IR-dark YSOs have 
luminosities that would impact on our analysis, we neglect 
such objects in this current work. 

2.3.3 Lyman flux 

In the RMS survey, massive YSOs are distinguished from 
young stars which have recently arrived on the MS by as- 
suming that the latter class of objects are ionizing their sur- 
roundings and driving a Hll-region. That is, sources with 
detectable Hll-regions are at a more advanced evolutionary 
state than YSOs. Therefore, in our model we must deter- 
mine the ionizing power of a star that is identified as having 
arrived on the MS. We then calculate the size and radio- 
brightness of the Hll-region, and then whether or not the 
Hll-region would be detected by the RMS survey. If the 
source has radio surface brightness above ImJy/beam at 
5GHz, then we assume that it wo uld be detected by the 
RMS continuum radio observations (|Urauhart et al. I l2007al . 




1.0 2.0 
Time (x10* yrs) 



Figure 4. The evolution of Stromgren radius Rss of a Hn-region 
around an accreting protostar, for two of the four different accre- 
tion scenarios studied here (with Mf{^=50MQ). The dotted lines 
illustrate the pressure driven expansion of the Hn-region once 
the star hits Mms(=2OM0), under the false assumption that 
the ionizing flux remains constant thereafter. The dashed lines 
show the growth of ijgs iii two different regimes: firstly, expan- 
sion at constant density due to the increase in ionizing photons 
as the star gains mass; and secondly, pressure-driven expansion 
once the mass has reached Mgn- The two approximations give 
different Hll-region sizes, and therefore give different radio fluxes. 
We make the simplifying assumption that the Hll-region size at 
any time is the larger of these two descriptions, illustrated by the 
solid line (offset for clarity). 



l2009bl ). If the Hll- region has a size larger than 2" in di- 
ameter, the source would be identified as being spatially 
extended in the high- resolution mid-infrared RMS imaging 
(|Mottram et al.ll2007l ). If the source does not exceed these 
limits, it is classified as a YSO. 

In calculating these quantities, we begin with a star with 
mass Mcur and determine the ze ro-age MS (ZAMS) stellar 
temperature and luminosity from lMevnet fc Maedei (|200G|) 
at this mas s. We then use the calc ulations of iMartins et al.l 
(2005) and lLanz fc Hubenvl (|2007l ). which provide the rate 
per unit area at which photons with frequencies above the 
Lyman limit Qlyc are emitted for a star of a given tempera- 
ture for O and B stars respectively. The relationship between 
stellar parameters and Lyman flux that we obtain is given 
in Table [T] 

There are several important points that should be noted 
re garding the numbers in T able [J Firstly, the cali brations 
of iMartins et al.1 l|2005l ') and iLanz fc Hubenvl (120071 ') do not 
provide us with Lyman fluxes for stars with masses above 
~6OM0 . Since the relationship between Lyman flux per unit 
area (<ZLyc) and stellar mass flattens off at high masses, we 
make the assumption that gLyc is constant above 6OM0. 
This will mean that we are likely to be slightly underesti- 
mating the Lyman flux of very massive stars. However, this 
will only impact the number of Hll-regions at high luminosi- 
ties (^ lO^ '^Z/0), which are very few in number in the -RMS' 
survey. 

Secondly, the absolute Lyman photon rates per unit 
time Qhyc as a function of stellar mass that we com- 
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Table 1. Table of ZAMS stellar parameters as a function of mass. The relationship between stellar mass, temperature, luminosity and 
radius is taken from.Mcync t fc Ma cdcr (2000). T he photometric p r opert ies BCv a-^d and the Lyman photon flux Qljc as a function 
of Teff are taken from lMartins fc Pled ||2006D and lLanz fc Hubenvl ll2007t) . and are corrected for differences in L*. t These values are likely 
to be slightly underestimated, see text for details. 



Mass 


Tcff 


Log (L*/Lo) 


R* 


BCv 


Mv 


LogQLyc 


tKH 




(K) 




(Rq) 






(phot s^i) 


(x lO^yrs) 


6 


19000 


3.01 


3.0 


-1.83 


-0.95 


43.33 


367.1 


9 


22895 


3.57 


3.9 


-2.28 


-1.89 


44.76 


175.5 


12 


26743 


3.96 


4.5 


-2.64 


-2.52 


46.02 


107.8 


15 


29783 


4.26 


5.1 


-2.88 


-3.01 


47.03 


76.2 


20 


33824 


4.61 


6.0 


-3.22 


-3.56 


48.00 


51.8 


25 


36831 


4.86 


6.7 


-3.47 


-3.93 


48.46 


40.2 


30 


38670 


5.02 


7.3 


-3.61 


-4.19 


48.69 


36.7 


35 


40596 


5.18 


8.0 


-3.75 


-4.45 


48.90 


31.7 


40 


42610 


5.34 


8.7 


-3.89 


-4.71 


49.09 


26.3 


50 


44589 


5.52 


9.8 


-4.03 


-5.02 


49.31 


24.1 


60 


46647 


5.70 


11.0 


-4.18 


-5.33 


49.43 


20.4 


70 


47662 


5.81 


12.0 


-4.25 


-5.53 


49. 51'!' 


19.7 


80 


48694 


5.92 


13.1 


-4.32 


-5.74 


49. 58'!' 


18.2 


90 


49449 


6.02 


14.1 


-4.37 


-5.92 


49. 65+ 


17.4 


100 


49913 


6.09 


15.0 


-4.41 


-6.06 


49.70t 


17.0 


110 


50381 


6.16 


16.0 


-4.44 


-6.21 


49.76t 


16.4 


120 


50853 


6.23 


17.1 


-4.47 


-6.36 


49.8lt 


15.5 



pute are different to those in iMartins et al.l ()2005l ) and 
iLanz fc Hubenvl (|2007l ). This is because these authors use 
a different caHbration between stellar mass and temper- 
ature to this study. Here, we use the relationship be- 
tween Mi,, Lj, TeS from the hydrostatic stellar models of 
iMevnet fc Maedej (|2000l ). 

Finally, we state explicitly that the final calibration of 
Qlyc versus M* should be considered to be uncertain at 
the 30-50% level, judging by the variation of numbers in 
the literature. However, the effect of this uncertainty on our 
results does not affect our conclusions, as will be shown later. 



H"*" to all levels with n >2. We then assume that one of two 
processes occurs: either the Hll-region undergoes pressure- 
driven expansion; or that the Hll-region grows in size at 
constant density due to the increase in Lyman flux from the 
accreting star, followed by pressure-driven expansion once 
the star reaches Mfln. We then determine which of these re- 
sults in the largest Hll-region. The differences between these 
two processes, in terms of the size evolution of the Hll-region, 
is illustrated in Fig. [l] 

For pressure driven expansion of an Hll-region (see e.g. 
iDvson fc Wilhamslll997l ). 



2.3.4 Size and radio emission of Hll-regions 

Accurately calculating the radio emission from a young star 
is a complex problem as there are many effects that one 
needs to consider: the Lyman fiux from the star, the den- 
sity and morphology of the ionized gas, the hydrodynamical 
expansion of the gas, and the expansion of the ionization 
region due to the increasing Lyman fiux from the accreting 
star as it grows in mass. However, accurate radio fiuxes and 
expansion rates are not the goal of this project, more that 
we want to get an estimate of the radio flux and the Hll- 
region size in order to determine how the source would be 
classified in the RMS survey. As such, we are able to make 
a number of simplifying assumptions when calculating the 
Hll-region preperties. 

Firstly, when the star's ionizing fiux is switched on, ei- 
ther at the end of the 'swollen star' phase, or as the mass 
of the star reaches Mms, an Hll-region is formed with an 
initial size given by that of a Stromgren sphere, 

"'-(IS?;)* <«' 

where rie is the initial electron density and /3b = 2.59 x 
10~^^cm"^s~^ is the radiative recombination coefficient of 




with Rs,o the Stromgren radius at t = tform; tcross = 
J?s,o/i'sound the sound crossing time of the Stromgren sphere 
at t = tform; ricfi the electron density at t = fform (assumed 
to be ~ lO'^cm"^), and Usound the sound speed in the gas 
(assumed to be ~10km s~^). 

We then determine the radio fiux S at a given frequency 
v of the Hll-region, 

5*. = B4T)[l-exp(-r,)] (19) 

where Bi,{T) is the Planck functio n, we assume a typical Hl l- 
region temperature of 8000K (e.g. lDvson fc Wilhamslll997l '), 
and Tv is the optical depth at frequency v, 

( V" X f^iZ^V" (20) 
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Table 2. Values of the parameters used in the fiducial model. 



Parameter 



Value 



YSO 



Yes 








Calculate size and 


radio flux of 


HII region 








Figure 5. Flow-chart for source classification. 



2.3.5 Source classification 

Once the observable quantities of each source in the simula- 
tion have been computed, we first apply the detection limits 
of the RMS survey to determine whether or not they would 
form part of the RMS sample. As the initial source selection 
was done using MSX point-source fluxes, the first detection 
limit to apply is that of the MSX Galactic Plane survey at 
21pm. As the completeness limits of the MSX survey are 
not particularly well defined, we ran our own completeness 
tests using sample images from the survey data to calcu- 
late the detection limits as a function of Galactic position. 
These completeness tests and their results are described in 
more detail in Appendix|B] For each source we compare the 
21pm luminosity to the completeness as a function of 21/im 
flux at the source's location on the sky, to determine the 
probability P that the source would be detected. A random 
number X is then generated, and the source is deemed to 
be detectable if X > P. 

Detected sources are classifled as either YSOs or Hn- 
regions based on their calculated continuum radio fluxes and 
angular sizes. Again, we must then compare these fluxes to 
the detection limits of the multiple radio surveys that were 
used in the compilation of the RMS survey. If a source is 
determined to have a radio surface brightness in excess of 
the 3(7-limit along that line-of-sight, it is classifled as a Hn- 
region. If the size of an object's Hn-region is greater than 
the spatial resolution of the imaging part of the RMS sur- 
vey (2"), it is assumed that it would be resolved and there- 



Accretion law 
rto 

Gas distribution 
IMF 

SFRcal 



McKee fc Tan (20031 
1.0 gcm~^ 
10*cm-3 
Spiral arms only 
Kroupa (2001) 
3M0yr~i 



fore again classifled as a Hn-region. If it is greater than the 
MSX beam size (~20"), it is assumed that the object would 
not belong to the MSX point-source catalogue, and so is 
discarded from the simulation. Objects which are not de- 
tected in the radio, are point-sources in the high-resolution 
RMS imaging, and are deemed to be detectable in the MSX 
survey, are classifled as YSOs. The classification process is 
illustrated by the flow-chart in Fig. [S] 



3 RESULTS 

In the following sections, we will describe the output of the 
model with various input parameters. Firstly, we concen- 
trate on the results when the 'flducial' set of input values 
are used. 



3.1 Fiducial model 

To flrst test the output of the model we start with a simula- 
tion for which the main parameters are set to fiducial values. 
These values are listed in Tabled 

For this initial model we start with the accretion law of 
MT03, with the parameter Eci set to 1.0 gcm"'^. In accord 
with the results of MT03, this implies that the mass at which 
massive stars join the MS, Mms, is set to 20Mq. The initial 
electron density rie when the Hn region flrst turns on is 
set to 10^c m~^, which is a typical val ue for ultra-compact 
Hn-regions l|Wood fc ChurchwellllT989l ). We use a Galactic 
gas distribution model in which star formation is conflned 
to the spiral arms only (i.e. the 'thin-disk' and 'thick-disk' 
diffuse gas components are ignored). This is because, in the 
NE2001 model, the thin-disk component includes a dense 
circumnuclear ring of gas at a Galacto-centric distance of 
i?GC~5kpc. This produces a large excess of sources at this 
location, which is not observed. The spiral arms on their 
own however reproduce the spatial distribution of sources 
reasonably well (see Sect. I3.2!2]) . 



3.1.1 niustration of classification criteria 

In Fig. [S] we illustrate the effect of our selection criteria on 
the output of the model. The top-left panel shows the 21/xm 
flux F21 of each object, demonstrating that no sources below 
the MSX detection limit (~4-10 Jy depending on location in 
the plane) are included in catalogue of young stars. The 
width and slope of the distribution of F21 as a function of 
Lboi are governed by the empirically derived parameters of 
Eq. (|14p . Another feature of the simulation illustrated by 
this plot is the absence of YSOs above Lboi ~ 10^ Lq. This 
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Figure 6. Various physical properties of the stars in a simulation using the fiducial parameters. In each panel the symbols have the 
foUowiing meanings: black-dot: all objects in the simulation; cyan-circles: objects classified as Hll regions based on their radio flux; 
blue-circles: radio-quiet objects classified as Hll regions based on their physical sizes; red-circles: YSOs; and yellow-circles: objects whose 
apparent sizes are larger than the MSX beam and so would be excluded from the RMS survey. The four panels show, clockwise from 
top-left: 21^m fiux as a function of luminosity; radio surface brightness as a function of luminosity; apparent angular size as a function 
of luminosity; and age of each object as a function of its current mass. 



corresponds to the ZAMS luminosity of a star with M* = 
AfMs(=2OM0 in this simulation). That is, any star which 
has reached the MS and is ionizing its nebula is found to be 
detectable by the RMS radio observations. 

In the top-right panel of Fig. |6] we plot the radio sur- 
face brightness of the objects in the simulation. The fig- 
ure shows that only sources with surface brightnesses above 
1 mjy/beam are detectable in the radio, and hence are classi- 
fied as Hll regions from their radio emission. Below this limit 
there are further sources which are classified as Hll regions 
based on their angular size: an object which is radio-quiet, 
but that would be found to be extended in follow-up infrared 
imaging with spatial resolution ~2" (i.e. has radius greater 
than 1"), is also classified as a Hll region. This is shown 
again in the bottom-left panel of Fig. [6l we also illustrate 
the point at which evolved Hll regions are discarded from the 
simulation, i.e. when their angular sizes become larger than 
18" and would no longer be point-sources in MSX. These 
'discarded' stars are shown as yellow points. 

Another aspect illustrated well by the top-right and 
bottom-left panels is that the observational and physical 
classification criteria match each other extremely well. Mas- 



sive YSOs are defined physically as proto-stellar objects 
which have not yet contracted to their MS configuration, 
emit very little Lyman flux, and so are not surrounded by 
ionized gas. Meanwhile, the observational definition is of a 
proto-star which is bright in the mid-IR but that shows no 
evidence of a Hll-region, either from continuum radio emis- 
sion of from spatially-extended mid-IR emission. In the top- 
right and bottom-left panels of Fig. [6] we see that there are 
very few objects misclassified as YSOs because their Hll- 
regions are too small or weak to be detected. 

3.1.2 Lifetimes 

The current age of the objects in each phase as a function 
of current mass is shown in the bottom right panel of Fig. 
[6] The objects classified as YSOs in the simulation all have 
masses below 20Mq (=Mms in this simulation). That is, 
as soon as a massive star reaches the main-sequence and 
begins to ionize its surroundings, it becomes immediately 
detectable in the -RMS' radio observations and hence is clas- 
sified as a Hll-region. 

The age distribution of YSOs is a function of cur- 
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Figure 7. The ages of objects in tlie simulation as a function 
of their final mass Mf{^. Symbols are the same as in the bottom 
right panel of Fig. |6] 

rent stellar mass. For objects in the simulation with 
Mcur<20MQthe YSO ages span l-3xl0^yrs. Some of these 
objects are those with high final masses which are actively 
accreting. However, most of these YSOs are objects which 
have finished accreting are still contracting to the MS from 
their 'swollen' phase (shown as green points). 

Once arriving on the MS, most objects are then picked 
up as Hll-regions. That is, there are very few YSOs outside 
the 'YSO envelope' clearly defined in the bottom right panel 
of Fig. E] by the green and red points (see also discussion in 
previous Section). Objects with masses below <1OM0 have 
very faint radio emission, but are still detected as Hll-regions 
as the angular size of the gas bubble they drive is larger than 
the spatial resolution of the RMS follow-up mid-IR imaging. 
The Hll-region phase itself spans around lO^yrs for lower 
mass objects, while it is about an order of magnitude shorter 
for the most massive stars. 

Another illustration of the phase lifetimes is plotted 
in Fig. [7] This plot demonstrates the time spent in each 
phase as a function of the final mass, Man. As stars initially 
grow in mass very slowly in the MT03 accretion rate model, 
no objects in the simulation become visible until around 
40,000yrs. The YSO phase then lasts ~ 1-3 x lOVrs for the 
lowest mass objects, a substantial fraction of which is spent 
contracting from the 'swollen' KH phase. For the most mas- 
sive objects, which can reach the MS while still accreting, 
the YSO phase lasts ~ 5 x lO^yrs. Therefore, in these simula- 
tions, though the most massive YSOs are currently ~2OM0, 
substantial numbers of these objects are still accreting and 
will go on to form stars well in excess of 20Mq. 

3.2 Comparison to the RMS results 

In the following sections we compare the results of our sim- 
ulations to those of the RMS survey. The features that are 
readily comparable are the relative distributions of objects 
as a function of luminosity and Galactic position. 

Within the RMS survey over 98% of the objects clas- 
sified as being YSOs or Hll-regions now have unambiguous 
kinematic distances. To assess the impact of the remaining 
2%, we compared three different observed luminosity distri- 



Figure 8. Relative numbers of YSOs and Ho regions as a function 
of stellar luminosity. The points with error bars show the results 
of the RMS survey, while the histograms show the results of the 
simulation with fiducial parameters. The thick red solid line and 
the thin blue solid line show YSOs and Hll regions respectively 
for the simulation with the fiducial parameters. The dotted lines 
show the results of a simulation in which the Kelvin-Helmholz 
contraction phase has been switched off. The real and simulated 
data have been binned identically, while the luminosities of the 
YSO and Hn observations have been offset from one another for 
clarity 

butions - one with the ambiguous distances set to the near- 
side distance, to the farside distance, and randomized be- 
tween near and far. In practice we found that the differences 
between these three distributions were within the noise, and 
so have a negligible impact on our conclusions. When com- 
paring with our simulated data, we adopted distances for 
these objects which were randomized between near and far. 

3.2.1 Luminosity distribution 

A powerful diagnostic of the RMS database is the relative 
numbers of YSOs and Hll-regions as a function of luminosity. 
The observational data are shown in Fig. [8] as data-points 
with error bars. It is clear that the observed distribution of 
YSOs peaks at a lower luminosity than does the Hll-region 
distribution, with m any more Hll regions found at higher 
luminosities (see also lMottram et al.ll201l"bh . 

Overplotted in Fig. [8] are the results of our simula- 
tion with fidicial parameters. We first note that our fidu- 
cial model already provides qualitatively good fits to the 
data. The quantitative fit to the Hll-region data is also very 
good, while the simulated YSO distribution predicts num- 
bers which are typically a factor of ~2 too large. As will 
be shown later, by fine-tuning the model parameters we are 
able to provide excellent fits to the data. 

The relative numbers of YSOs and Hll-regions (i.e. 
iV(YSO)/Af(Hll)) are also well matched. This is a valida- 
tion of the approximations we made in calculating the ob- 
served properties of the Hll-regions. Under our current as- 
sumptions, almost every object in the simulation which has 
evolved as far as the MS is detected as a Hii-region (see also 
discussion in Sect. l3.l"T|l . If our assumptions underestimated 
the radio flux Si, then the numbers YSOs/Hll-regions would 
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be unchanged, as nearly all MS objects are already detected. 
Meanwhile, if we overestimated Si^ then the Hll-regions 
would still be picked up as being spatially extended (see top- 
right panel of Fig. [6]). The physical sizes of the Hll- regions 
are goverened by their expansion rates, which in turn is dic- 
tated by the parameter rio. This can greatly alter the num- 
bers of Hll-regions in the simulation while leaving the num- 
bers of YSOs unchanged. Using the typical ultra-compact 
Hll-region value of nc=10''cm~^ (|Wood fc Churchwell|[l989l ') 
we find excellent agreement between the simulated luminos- 
ity distributions of YSOs/Hll-regions and the RMS results. 

One obvious feature of the simulated data is the steep 
drop-off in YSOs at log(L/LQ) = 4.9, whereas the obser- 
vations show a smoother drop-off at high luminosities. The 
location of this drop-off corresponds to the luminosity of ob- 
jects with large Mfln when they arrive on the MS at Mms- 
The parameter Mms depends on the accretion rate - higher 
accretion rates mean that stars join the MS at higher masses. 
Therefore, we can adjust the location of the high-L* dropout 
by changing Mms, or equivalently the accretion rate. 

Finally, we address the topic of the pre-MS 'swollen' 
phase of massive stars. In Fig. [S] the solid and dotted lines 
show the luminosity distributions of simulations with and 
without a KH contraction phase respectively. As we have 
already illustrated in Figs. [6] and [7] the inclusion of a KH 
phase significantly prolongs the YSO phase in lower mass 
objects, and results in larger numbers of low-i* YSOs. In 
Fig. [S] we see that if we disregard the KH phase the num- 
ber of low-L* YSOs is reduced. Meanwhile, the number of 
high-L* YSOs remains approximately the same, since any 
KH phase for these objects would be short compared to the 
age of the star when accretion terminates. The numbers of 
Hll-regions remain unchanged, since in our simulation the 
duration of the Hll phase is the same whether or not it is 
proceeded by a KH phase. These effects result in a poorer 
qualitative fit to the data compared to a simulation which in- 
cludes a KH phase. We therefore conclude that a KH swollen 
phase is required for stars with masses <2OM0in order for 
our simulations to provide a good match to the data. 

3.2.2 Galactic distribution 

Another diagnostic of the RMS sample is the distribution of 
sources with Galactic position. As well as comparing with 
the source distribution in Galactic longitude I and latitude 
b, the radial velocity of each RMS source gives us an unam- 
biguous measurement of Galacto-centric distance -Rgc, as- 
suming that the sources follow the Galactic rotation curve 
to within a few km s~^. 

We again state that in our fiducial model we consider 
the star formation to follow the gas distribution in the spiral 
arms only. When including the full gas distribution model, 
which has components due to the thick disk as well as a ring 
of gas in the thin disk at _RGC~5kpc, we found an extremely 
large overdensity of sources at low Rac which is not seen in 
the data. Our results suggest that massive star formation is 
confined to the spiral arm component of the Galaxy. 

The total distribution of YSOs and young Hll regions 
from the RMS survey with Galactic position and compar- 
isons with the model results are shown in Fig. [9] As with the 
luminosity distributions in the previous section, the approx- 
imate normalization and overall trends of the simulations 



match the data rather well. However, there are qualitative 
features of the simulation which are not matched by the 
observed data. Firstly, there are two 'horns' in the number 
counts of the model's /-distribution, at I ~ 50° and I ~ —70°. 
Secondly, in the 6-distribution the fall-off in source counts as 
one moves away from the plane is steeper in the observations 
than in the simulation. Thirdly, the right-hand panel of Fig. 
[9] shows that there is a surplus of sources in the simulation 
at _RGC~8kpc with respect to the data. 

These three features can all be understood as being due 
to the same feature in the Galactic gas distribution model. 
The horns in the Z-distribution are due to an overdensity of 
sources in the Sagittarius-Carina arm. As this arm passes 
close to the Sun, this also causes the overdensities of sources 
at high Galactic latitudes and at i?GC~8kpc. In principle 
the RMS results could be used to fine-tune the description 
of the Galaxy's gas content, however this is not the goal 
of the current work and we leave such a study for a future 
paper. 

3.3 Effect of varying fiducial parameters 

3.3.1 Stellar IMF and SFRg^i 

Broadly speaking, these parameters control the overall nor- 
malization of the predicted number of objects in the RMS 
survey. The IMF of Krouoa (2001) puts a significant amount 
of mass into objects with sub-Solar masses. If this were to 
be changed to another description of the IMF, such as a 
standard Salpeter law down to O.SMq, this would increase 
the number of massive objects in the simulation, though the 
relative numbers of objects as a function of mass would be 
unchanged. Though there is speculation as to variations in 
the slope of the IMF at high masses and t o the presence of 
an upper-mass cutoff (jBastian et al.l[2010l ). the numbers of 
objects in the RMS survey with inferred high masses (i.e. 
those with Lboi > IO^Lq) are too few for us to test this. 

The parameter SFRcai linearly affects the overall nor- 
malization, for example doubling SFRcai results in twice 
as many objects in the simulation. Changing the IMF to 
a Salpeter IMF would result in more massive stars, which 
would then require SFRq^i to be decreased to maintain the 
same overall number counts of objects. However, the major- 
ity of independent estimates of SFRcai are above 2MQyr~^, 
so there would be little justification in reducing this param- 
eter beyond this value. The consistency in overall numbers 
between our simulations and the observations serve to jus- 
tify our initial choices in SFRd^i and the shape of the IMF. 
However, later in this study we will allow some variation in 
SFRos^i in order to fine-tune the overall normalization. 

3.3.2 Initial Hll region density Ue 

As has already been discussed, changes to this parameter 
affect only the duration of the Hll region phase. Increas- 
ing rio by a factor of ten from the fiducial value (i.e. to 
rio^lO^cm"^) causes the Hll regions to be initially more 
compact, and expand at a slower rate. As a result, the num- 
ber of Hll regions in the this model is increased by a factor 
of ~5 above those of the fiducial model. Similarly, dropping 
the density to no=10'^cm~^ causes the Hll regions to rapidly 
expand once switched on, shortening the phase dramatically. 
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Figure 9. The distribution of YSOs and young Hll regions with Galactic position. The blue data-points with error bars show the results 
of the RMS survey, while the histograms show the results of the simulation using the fiducial model. 



leading to an under-prediction of the number of Hll regions 
again by a factor of ~5. The fiducial value of ne=10^cm~'^ 
gives a good fit to the observed Hll-region distribution across 
all luminosities. 

We note that, in reality, Hll-regions are likely to be 
expanding into a density gradient, leading to cometary mor- 
pholgies wi th size evolutions that differ from those calcu- 
lated here ([Arthur fc Hoard l2006l ). Our use of a constant 
He simply represents an average over the early stages of a 
Hll-region's evolution. 

3.3.3 Stellar mass - Qhyc relation 

As discussed in Sect. 12.331 the relationship between stellar 
mass and the i onizing flux Qlvc is u ncertain. As has already 
been noted bv lMartins et al] (|2005l ) the uncertainty in Qlyc 
for a star of a given mass could be in excess of 50%, though 
we find from comparisons of various authors' estimates of 
the mass - spectral-type relation that the uncertainty in 
Qhyc for a given spectral-type may be as high as ~0.7dex, 
or roughly a factor of five. 

To investigate the effect of this uncertainty, two further 
simulations were run with Qlyc scaled by ±0.7dex. With 
Qhyc decreased, Hll regions emit less radio fiux and expand 
at a slower rate, meaning that the Hll region phase is pro- 
longed. However, the vast majority of Hll regions are still 
picked up by the selection criteria, meaning that very few 
Hll regions are classified as YSOs in the simulation. The pri- 
mary effect of decreasing Qhyc is therefore to increase the 
number of Hll regions, though only by a factor of ~2 in most 
luminosity bins. 

Similarly, increasing Qhyc by 0.7dex has only a mild 
impact on the results. The overall number of Hll regions is 
lowered due to the shorter Hll region lifetime, though again 
only by a factor of ~2. So, despite the large uncertainty in 
the Qhyc-Mcur relation, we conclude that this uncertainty 
has negligible impact on the results of our simulations. 

4 DISCUSSION: A COMPARISON OF 
DIFFERENT ACCRETION LAWS 

Thus far, all our models have used the MT03 prescription of 
accretion, in which AfcurOC M^^t^, i.e. accretion which ac- 



celerates with time. In the following sections we investigate 
the free parameters in this prescription, as well as test other 
modes of accretion. 

We measure how well a certain accretion law fits the 
data from the YSO Lboi-distribution only. This is because 
of (a) the large number of assumptions that have gone in to 
predicting the Hll-region properties, and (b) the morphol- 
ogy and normalization of the Hll-region Lboi-distribution 
can be augmented by altering parameters such as rio and 
Qhyc, which are uncertain to factors of ~five. For now, we 
consider it a success of our model that the Hll-region Lboi- 
distribution is matched very well for the fiducial parameter 
values. 

To assess the results produced by different accretion 
laws we employ a standard reduced-x^ analysis between the 
observed and simulated Lboi-distributions. We have adap- 
tively rebinned the observed YSO distribution such that 
there is a minimum of 9 objects per bin, i.e. that each bin 
has a minimum fractional error of 33% and therefore has a 
minimum significance of 3a in Poissonian statistics. For the 
simulated data we re-ran each simulation 50 times to reduce 
random errors to well below those of the observations. 

We allowed the precise normalization of the Lboi- 
distribution to vary, since the parameter that affects overall 
numbers of objects (i.e. S-Fi^cai) is itself uncertain, and we 
would not wish to discriminate against a particular accretion 
law simply because of a small offset in numbers. However, 
where a large normalization factor was required to achieve 
a minimal we critically assess the implications for the 
associated re-scaling of SFRqi^i. The best-fitting a-nd as- 
sociated SFRqs^i values for each model run are summarized 
in Table O 



4.1 MT03 accretion 

In the MT03 turbulent core model, the accretion rate of a 
star of a given final mass is dictated by the clump pressure 
and pressure structure. In their fiducial model this is pa- 
rameterized in terms of the clump's surface density Ed. In- 
creasing the accretion rate not only leads to faster formation 
times, but also results in high Mgn stars joining the MS at 
higher masses (i.e. Mms increases, due to greater 'swelling' 
in the accretion phase). According to MT03, varying Ed be- 
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Table 3. Summary of each accretion scenario and associated 
model parameter, its value, and the implied average star- 
formation rate of the Galaxy. 
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Figure 10. The effect of varying Sd, or equivalently the time- 
avcraged accretion rate a function of final stellar mass, on the 
simulated luminosity distribution of YSOs. The histograms show 
the simulation results for values of S(.i=0.1, 1.0 and lOgcm"^, 
which have been renormalized to best match the data by scaling 
the overall star- formation rate (see text for details). The data 
have been adaptively binned to have a minimum error per bin of 
33% (i.e. a 3a detection in Poissonian statistics). 



tween 0.1 and lOgcm" results in values of Mms between 
12 and 24M0 . In Fig. [10] we show luminosity distributions 
using the extrema for these values of Ed and Mms- 

The low-Eci model has large differences to the fidu- 
cial model. The decreased accretion rates have the effect 
of lengthening the YSO phase, which is now dominated by 
stars still accreting rather than by stars which have finished 
accreting and are contracting to the MS. This means that 



the total number of low luminosity YSOs in the simulation is 
much larger, and to match the observed YSO levels in Fig. llOl 
we had to normalize the data by reducing the numbers of ob- 
jects by a factor of 6. The only physical explanation for this 
normalization factor would be to reduce the Galactic star 
formation rate to S'7<'7?Gai~O.5M0yr~^, whereas the vast 
majority of estim ates of this quantity are >2M0yr~^(see 
iDiehl et al.l |2006| , and refs therein). Also, the decrease in 
Mms pushes the YSO drop-out to lower luminosities, mak- 
ing it much more pronounced, and in clear disagreement 
with the RMS results. This is reflected in the poorer 
value when compared to other model runs of this scenario. 

The fiducial and high-Eci models match the observa- 
tions very well. The high-Ed model has a higher value of 
Mms than the fidicial model, the effect of which is to push 
the YSO drop-out to higher luminosities. As the steepness of 
the IMF dictates that there are fewer objects at these higher 
luminosities, this makes the drop-out slightly less dramatic. 
The faster accretion rates mean that massive stars form 
quicker, however this does not greatly affect the distribu- 
tions of YSOs. The YSO population is dominated by those 
objects with A/fln<MMS, and a substantial fraction of the 
YSO lifetime of these objects is the Kelvin-Helmholz con- 
traction phase, which is independent of the accretion rate in 
our model. The objects with Mfln>MMS spend very little 
time in the YSO phase in the fiducial model (cf. Fig. [7]), 
so decreasing the contribution of such objects to the YSO 
population does little to the observed YSO distribution, as 
the high Mun stars spend most of their formative years as 
Hll regions. 

In terms of which value of Ed produces the best fit, 
the fiducial model achieves the lowest reduced x^- the 
Ed=lgcm^^ model has x^ = l-2, and the Ed=10gcm~'^ has 
X^=1.4. Given the similarity of these values, and the random 
uncertainty associated with how the data are binned, we 
consider each of these models to produce equally good fits. 
However, the fiducial value of Ed is comparable to that typi- 
cally observed in Galactic star-forming clumps l|Plume et al.l 
1 19971 ). whereas a value of Ed = 10gcm ^ seems excessively 
high and is not supported by observations. The renormal- 
ization factor for the Ed^lgcm"'^ model is 0.65, implying 
a global star formation rate of S'f7?Gai = l-6M0yr~^. This 
is reasonably consistent with the vast majority of estimates 
of this parameter, which place it in the range 2-4M0yr~^. 
Therefore, we consider that the Ed = lgcm~'^ model pro- 
vides a good match to the data while being consistent with 
complimentary observations of star-forming clumps. 



4.2 Bondi-Hoyle accretion 

The simple Bondi-Hoyle (BH) 'competitive' accretion model 
we are using has the attractive property that it has very 
few free parameters o nce we employ the normalization of 
iBonneU fc Bate! (|2006D . We are not aware of any robust nu- 
merical simulations of BH accretion which predict properties 
such as Mms and MS contraction time, so in these simula- 
tions we have assumed that the MS contraction time ap- 
proximates to the Kelvin-Helmholz time, in accord with the 
results of HO09. We allow Mms to remain a free parame- 
ter, and explore the range Mms=15-3OM0 which are typical 
values taken from the other accretion laws studied here. In 
Fig. [TT] we show the results of these simulations with three 
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Figure 11. Same as Fig. llOl but for Bondi-Hoyle accretion mod- 
els, for three different values of the mass at which accreting stars 
arrive on the main-sequence, M^s- AH simulations have been 
re-normalized by a factor of 1.5, see text for details. 



Figure 12. Same as llOl but for uniform constant accretion rates. 
Each histogram has been re-normalized to fit the observed num- 
bers of low-Lboi YSOs, see text for details. 



diflterent values of Mms between 15 and 30Mq. The qualita- 
tive agreement with the RMS data is very good when Mms 
is between 20-30Mq, while the overall normalization is also 
very good: optimal renormalization factors were found of 
-0.6, implying SFRq^i between 1.7-2.OM0yr"^ 

The best fitting model is that with Mms=20Mq, with 
reduced x^=l-0. This value of Mms seems reasonable, 
given the quantitative results of MT03 and HO09 who find 
Mms = 12-3OM0 depending on the accretion history of the 
protostar. The optimal normalization of this model requires 
SFRci,i=1.9MQyr~^ , which again is reasonable given other 
measurements of this value. 



4.3 Constant accretion 

In this section we explore two distinct scenarios. Firstly, we 
investigate a uniform accretion rate M for all stars, regard- 
less of mass. Secondly, we fix the accretion rate such that 
Affin/M is constant, that is that all stars finish accreting on 
the same timescale tform, independent of Man. 
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Figure 13. Same as 110! but for uniform formation timescales 



and with accretion rates which are constant with time. 



The histograms show simulations with tfomi of 10 i 3 X 10 and 
7 X lO^yrs, and have been multiplied by factors of 0.4, 0.1 and 
0.2 respectively to fit the observed levels. 



4.3.1 Uniform M 

The luminosity distributions for three values of M are shown 
in Fig. 1121 using the computed tracks at Af =10~^MQyr~^ 
and 10~*A^Qyr~^, and a track linearly interpolated between 
these two at lO~^'^Af0yr~^. For each track we found the 
value of AfMS by plotting the track on a H-R diagram and 
finding the approximate mass at which the track joined 
the ZAMS of ( Mcynct & Macdcr 20o3). The value of Mms 
was Ri27M0 for the lO"^'^ and lO"''M0yr"^ tracks, and 
«36M0 for the lO"^M0yr"^ track. 

While the normalization is reasonable (see Table [3] for 
the implied values of S-F7?Gai), these models give poor qual- 
itative fits to the data. The reason for this is that in the 
HO09 accretion tracks, the behaviour of Lboi with Mcur is 
such that at some point Lboi reaches a temporary plateau. 
In the simulation, this causes a pile-up of objects at a cer- 
tain luminosity. For example, in the Af=lO~^M0yr~^ this 



plateau occurs at Lboi— 1O'''^L0, leading to a spike in the lu- 
minosity distribution at this value (see the blue dashed curve 
in Fig. [HI). Similarly, in the M=lO"^M0yr"^ model, this 
pile-up occurs at Lboi— 1O^'*L0. The interpolated model 
(Af=lO~'''^Af0yr~^) does not have any such pronounced 
spikes in the luminosity distribution as the interpolation 
process serves to smooth any such spikes out. However, this 
model still produces a poor qualitative fit to the data, the 
Lboi-distribution being significantly flatter than the obser- 
vations. 

4.3.2 Uniform tform 

Here, we experiment with using accretion rates which are 
constant with time, but depend on the final mass of the 
star in such a way that all stars accrete their matter on 
identical timescales. To do this, we took the birthlines of 
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HO09, which are sampled at accretion rates of 10~^, 10~* 
and lO~^M0yr~^ and linearly interpolated them onto a finer 
grid, sampling at every O.Oldex. For each star in the simula- 
tion we then calculated the accretion rate required for it to 
form in a time tform, where tform is a free parameter, and as- 
signed each star a constant accretion rate of M=Mfln/iform- 
The lowest values of M were clipped at lO~^Af0yr~^, and 
the highest at 10~''MQyr~^, as we have no birthline calcu- 
lations outside this range. 

The point at which each star joins the MS, Mms, is 
then a function of the star's final mass. To determine Mms 
as a function of M we plotted each interpolated track on a 
H-R diagram and identified the approximate point at which 
it joined the ZAMS track of (|Mevnet fc Maedeij[200oh . Stars 
with A'/fln<A/MS were again assumed to undergo a contrac- 
tion phase before arriving on the MS, with the length of this 
phase equal to the K-H timescale of a star of that mass. 

We can summarize the predicted behaviour of such 
models as follows. The effect of decreasing iform is to increase 
the accretion rates. This means that stars can pass through 
the YSO phase more quickly, resulting in less YSOs in each 
luminosity bin. However, this effect is mitigated by the fact 
that the higher accretion rates push the value of Mms to 
higher masses, extending the length of the YSO phase for 
objects with high Man- Altering tform can therefore affect 
both the qualitative shape and the overall normalization of 
the Lboi-distribution. 

The results of simulations with three different values of 
uniform tform are shown in Fig. 1131 As expected, in the sim- 
ulation with short formation timescales (tform=70 x lO'^yrs) 
the YSOs have higher luminosities than in the simulations 
with longer tform. The large value of A4ms (~4OAf0) means 
that the numbers of YSOs with Lboi<; IO^Lq are overpre- 
dicted by factors of 2-3 compared with the lower luminosity 
objects. At longer formation timescales of fform=3 x lO^yrs 
the YSO drop-out matches the observations slightly better, 
though the simulated distribution is again 'top-heavy' com- 
pared to the observations. The simulation with the longest 
iform produces a very good qualitative fit to the data. How- 
ever, the predicted numbers of YSOs are too large by a fac- 
tor of five, requiring SFRcai to be reduced to O.GAfQyr"^ to 
match the observed levels. 
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Figure 14. Same as llOl but with accretion rates from SK04 wliich 
decelerate with time. Simulations are shown using three different 
values of M ms • The simulations have been renormalized to fit the 
data by reducing the total number of objects in the simulation 
by a factor of three. 



els. In Fig. [14] we show simulation results for three different 
values of Afivis- 

The qualitative shape of the YSO Lboi-distribution this 
time does not match the observations. The longer time that 
YSOs spend at higher masses compared to the accelerating 
accretion models produces a luminosity distribution which 
is more top-heavy than the observations. This is reflected in 
the values, which are higher than those of the accelerating 
accretion models. 

The overall numbers of objects, as with the constant 
rate accretion law, are overpredicted in this model - to 
achieve the correct normalization each histogram had to be 
scaled down by factors of 2-3, implying global star-formation 
rates below l-1.5Af0yr~^. Meanwhile, the qualitative shape 
of the models' Lboi-distributions is much flatter than the 
observations. This is due to the very high accretion rates 
that objects have at early times, and hence higher accretion 
luminosities. The net effect is that lower mass objects are 
moved to higher luminosity bins with respect to accelerating 
accretion models, producing a flatter Lboi-distribution. 



4.4 Decelerating accretion 

As already discusse d in Sect. [H t he prescription of deceler- 
ating accretion from lSchmeia fc" Klcssen (2004) cannot pro- 
duce massive stars within a reasonable formation timescale 
when using the standard values for the constants in Eq. pip . 
Massive stars can be formed however if the constants are 
pushed to the boundaries of their uncertainties. This max- 
imises the value of r in Eq. H13p . and allows the protostars 
to spend more time at high accretion rates. For the analy- 
sis in this Section, we use values of tq — 4.9 x lO'^yrs and 
n = 6.3 X 10*yrs. 

Figure O shows the luminosity distribution of young 
massive stars from the simulation using the decelerating ac- 
cretion model. The SK04 accretion model provides no esti- 
mate of Mms, but the low accretion rates as the star is ap- 
proaching its flnal mass suggest that A^ms may be lower un- 
der this model than for the constant and accelerating mod- 



4.5 Which accretion laws fit the data best? 

We now discuss which of the accretion models we have tested 
best fits the observed luminosity distribution of the RMS 
survey data, taking into account the results described in 
this section thus far. 

The first conclusion that we draw is that the deceler- 
ating accretion scenario is the least likely of all those stud- 
ied here. While it may be argued that the prescription of 
SK04 is only calibrated between l-lOAf0, the poor match 
between the data and the observations is still evident at low 
luminosities (and hence at masses between 5 and 2OAf0). 
Also, irrespective of the prescription's quantitative details, 
the qualititave effect of a decelerating accretion law is to 
produce a flat luminosity distribution of YSOs. Therefore, 
we conclude that this mechanism is unlikely to be the preva- 
lent mode of massive star formation in our Galaxy. There 



Critical Tests of Accretion Models for Massive Stars 17 




Simulation: RIVISdata: 

HII X 

YSO A 



4.0 



4.5 



5.0 5.5 
log(L/L^ 



6.0 



6.5 



Figure 15. Same as Fig. |8]but with the best-fit model, using the 
BH accretion law with Mms=23M0. To match the Hll-regions, 
we used an initial electron density of nc=4 X lO^cm^. 



are however caveats to this statement as we made several 
approximations in our calculations, such as that of the accre- 
tion luminosity of massive protostars, and of objects joining 
the MS mid-way through their accretion phase just as the 
accelerating accretion models do. Stronger critical analysis 
of the decelerating accretion scenario awaits the numerical 
predictions of the accretion luminosities and pre-MS tracks 
of massive stars forming in this way. 

Similarly, we also argue against the uniform accretion 
rate models. Models where all stars accrete at exactly the 
same rate produce spikes in the luminosity distribution due 
to the discrete features in the HO09 accretion tracks, while 
the global gradients of the simulated Lboi-distributions are 
much flatter than that observed. 

The uniform fform models do produce very good quali- 
tative fits to the data, but only if the formation timescale is 
very long, i.e. ~ 7 x lO^yrs. As previously discussed, using 
the fiducial value of SFRq^i this produces greater numbers of 
YSOs than are observed. To recover the observed numbers, 
it requires that SFRqi^i be reduced to O.6M0yr~^, a factor 
of five lower than the typically-quot ed value of SMpyr"^, 
and b elow even the lower limit of Robitaill e fc Whitney! 
l|20ld ). whose estimate of SFRca.! is significantly lower than 
any other recent measurements. Further, this formation 
timescale would imply accretion-rates for massive stars of 
W^Mqyt-^, which are an order of magnitude lower 
than current observa t ional and theor etical estimates (e.g. 
iKrumholz et al. l l2009l : lQiu eralll201ll ). 

The best fits are produced by the two accelerating ac- 
cretion models. Both produce reduced statistics close 
to unity (see Table O, and hence provide excellent fits to 
the data given the experimental uncertainties. At the same 
time, the number of sources in each model is a close match 
to the observations, meaning that the implied SFiJcai re- 
mains sensible: we find best fitting average Galactic star- 
formation rates of l.BMoyr"^ and 2.OM0yr"^ for the MT03 
and BH models respectively. These are co nsistent with most 
recent measurements of this property fsee lDiehl et al.|[2006l : 
iRobitaille &: WhitnevlfioiOl . and references therein). 

At this stage, we are unable to conclusively state which 



of the accelerating models fits the data best. Random er- 
rors in, for example, how the data are binned, mean that 
the difference in the x^statistic between the two models 
is not significant. We should state that the MT03 model 
has a more rigorous quantitative basis compared to our 
BH model, which simply assumes the analytical formula for 
Bondi-Hoyle accretion and ad- hoc estimates for Mms • Thor- 
ough numerical calculations of each scenario may in future 
allow us to distinguish decisively between the two. 

In Fig. [15] we plot the YSO and Hll-region luminosity 
distributions for the best-fit model. We have chosen the BH 
model, since this gives a slightly lower than the MT03 
model, and we have increased Mms to 23M0 to match 
the largest YSO luminosity bin, though this is more for 
aesthetic purposes since the reduced is no better than 
for Mms=2OM0. Since the overall level of Hll-regions was 
slightly underpredicted, we have increased the initial elec- 
tron density to nc=4 x lO^cm"^ to slightly prolong the Hll- 
region phase and therefore increase the overall number of 
objects. Note that this value of ric is consistent with the 
typically quoted densities for ultra-compact Hll-regions of 
no> lO-'cm^ (Wood & ChurchwcU 1989). 

We can highlight two features of Fig. [15] where there 
is room for improvement in the fit. Firstly, the simulated 
YSO cutoff occurs at lower luminosities than in the observed 
trend. A better fit may be obtained by using values of Mms 
which depend on Mgn rather than a single blanket value. 
Since Mms depends on accretion history (see HO09, plus 
discussion in this work), which in turn depends on stellar 
mass in the accelerating accretion rate models, it is entirely 
plausible that Af ms may be dependent on final stellar mass. 
This would then serve to smooth-out the YSO cutoff and 
make it more like the observations. 

Secondly, there is an apparent downturn in the Hll- 
region luminosity distribution at low luminosities, which is 
not reproduced in our simulations. An improvement to the 
fit may be obtained if lower Mfln objects have lower initial 
gas densities once they have reached the MS. Since low den- 
sity regions expand at a faster rate, this will cause the Hll- 
region phase lifetime to be shorter for these objects. Though 
there are a large number of assumptions that go into our 
simulation of the Hll-region distribution, it is not an unrea- 
sonable expectation that higher mass stars should form in 
regions of higher gas densities. 

4.6 Phase lifetimes of the best-fit model 

For our best-fitting model described in the previous Section, 
we now revisit the topic of the YSO and Hll-region phase 
lifetimes. For this, we compare to our rece nt empirical es- 
timates presented in lMottram et al] l|2011al ). In that study, 
the lifetimes were estimated by dividing the number of ob- 
jects per luminosity bin by the number of OB stars with 
the same luminosity, then multiplying by the MS lifetime 
of OB stars. There are three main sources of uncertainty 
in estimating phase lifetimes in this way. Firstly, the OB 
star population of the Galaxy is only complete within a few 
kiloparsecs of the Sun, which serves to increase the random 
errors. Secondly, the relationship between spectral type, stel- 
lar luminosity and stellar mass is still uncertain, and may 
introduce systematic errors, though we mitigate this uncer- 
tainty by using the same sources for these relations. Finally, 
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Figure 16. The phase lifetimes, predicted by the best-fit model, 
for both YSOs and Hll-regions, as a function of Lbol- The dot- 
ted and dashed lines in dicate the lifetimes derived empirically by 
iMottram et alj bOllal) . 

this method gives phase lifetime as a function of luminosity, 
which for actively accreting objects is variable throughout 
their formation. However, our simulation provides a 'snap- 
shot' of the Galaxy's population of young massive stars in 
a way designed to mimic the RMS survey. Therefore, we 
can make a quantitative measurement of lifetime versus Lboi 
from our simulation, despite the fact that objects may mi- 
grate from bin to bin throughout their evolution. 

To calculate the lifetime as a function of luminosity, we 
first construct a diagnostic diagram of the simulation sim- 
ilar to that of the lower-left panel in Fig. [S] but plotting 
age against Lboi- Then, for a given luminosity bin, we deter- 
mine the age range (i.e. the difference between the minimum 
to maximum ages) of all YSOs/Hll-regions in that bin, dis- 
carding the oldest and youngest 5% of objects to limit the 
impact of statistical outliers. The results of this analysis are 
pl otted in Fig. | 16l where we also compare with the results 
of lMottram et al.l (|2011al ) which are illustrated by the long- 
and short-dashed lines. 

We find that the lifetimes of YSOs are typically 
lO^yrs with shorter lifetimes for the more luminous objects, 
predominantly an effect of the reduction in KH timescales 
at higher masses. On the other hand, Hll-region lifetimes 
show a weaker dependence on luminosity, spanning a factor 
of ~3 between across the plotted luminosity rangc^. This 
can easily be understood from a numerical analysis of Eqs. 
()16p and (|18|) and the Lyman fluxes from Table [1] (see Sect. 
I2.3.4|) . which show that the expansion rates of Hll-regions 
span a similar factor of ~3 across the range of luminosi- 
ties plotted in Fig. 1161 under the assumption of initial gas 
densities which are uniform across all objects. 

The agreement between the two estimates of Hll-region 
lifetimes is excellent. Recall that we have already fine-tuned 
the initial electron density Uc in order to match the observed 
numbers of Hll-regions. Since the observed number of objects 
at any given time is a product of the production rate and the 

Note that this lifetime refers to that of compact phase of the 
Hll-region, such that it would be detected in the RMS survey. 



lifetime, the fact that we can match the lifetimes and overall 
numbers simultaneously must mean that we also have an 
accurate estimate of the production rate (i.e. the globally 
averaged star- formation rate). In other words, the overall 
number of Hll-regions depends on SFTJcai and n^, whereas 
the lifetime depends on only ric. 

In terms of the YSOs, there is an apparent discrepancy 
be tween the two e stimates. The systematic errors quoted 
bv lMottram et al] (j2011a), ±0.45dex in the zero-point and 
±0.5 in the gradient, mean that the differences are signif- 
icant only at the l-2a level. In addition, the differences 
in gradient can in par t be explained by the way in which 
iMottram et al.l l|2011al 'l determined the YSO lifetimes. As 
discussed above, these authors assumed a 1:1 correspon- 
dance between the luminosities of the YSOs and the those 
of main-sequence (MS) stars. However, since YSOs are often 
accreting objects they will have some accretion luminosity, 
and so their total luminosities will exceed those of MS stars 
with the same mass. In contrast, the majority of Hll-regions 
in our simulation have finished accreting, and have bolomet- 
ric luminosities comparable to their MS luminosities. 

For the simulation used in Fig. 1161 we found that if the 
YSOs bolometric luminosities were used to estimate Mcur 
(i.e from interpolation of Table [T|, then the current masses 
were typically overestimated by around 30%. This leads to 
the ratio of the number of YSOs to the number of MS stars 
at a given mass being over-predicted, and ultimately to over- 
estimates of the YSO lifetimes. Quantitatively, we find that 
this overestimate in lifetimes is ~70% at Lboi = lO*I/0, and 
40% at Lboi=lO^Z/0. Though this effect is not large enough 
to fully explain the discrepancy in YSO lifetimes in Fig. 1161 
it may be a contributing factor to the difference between the 
two gradients. 



5 SUMMARY 

We have constructed a model to simulate the distributions 
and physical properties of massive protostars throughout the 
Galaxy, with the aim of predicting the observed luminos- 
ity distribution of massive Young Stellar Objects (YSOs) 
from the .RMS' survey. To compute the observed properties 
of each protostar in the simulation we employed several dif- 
ferent prescriptions of the rates at which stars accrete mass 
as a function of time. For each accretion model we identi- 
fied the parameters which produced the best qualitative fit 
to the observed YSO Lboi-distribution, allowing the global 
star formation rate SFT^cai to vary to optimize the overall 
normalization. 

Our main findings may be summarized as follows: 

• The luminosity distribution of YSOs is best described 
using accretion rates which increase as the star grows in 
mass. Accretion rates which were constant in time or which 
decreased as the star grew in mass were ruled out, as they 
either produced poor qualitative fits, produced far too many 
YSOs, or both. 

• The lack of YSOs at high luminosities is consistent with 
a picture of star-formation whereby stars arrive on the main- 
sequence (MS) once they reach a mass of 2O-3OM0, even if 
the final mass of the star is well in excess of this value. We 
suggest that the precise mass at which stars arrive on the 
MS may be some function of the final stellar mass. 
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• In order to produce a satisfactory qualitative fit to the 
YSO distribution at low luminosities, an intermediate radio- 
quiet phase is required between the end of accretion and the 
ignition of a Hll-region. The length of this phase is well de- 
scribed by the star's Kelvin-Helmholz timescale, and there- 
fore is consistent with a 'swollen star' period of pre-MS evo- 
lution for stars with masses ~5-2OM0. 

• The maximum final stellar mass Man at which stars 
no longer experience this 'swollen' phase is indicated by the 
observed drop-off in YSO numbers at Lboi ~ This 
corresponds to stars with masses Mfln^SOM©, which in our 
simulations have accretion rates of 1-4 xlO~*M0yr^^. 

• From the best-fit models we find that the overall num- 
bers of objects can be well reproduced by a globally averaged 
star-formation rate of the Galaxy of 1.5-2M0yr~^. 

• Our best-fit model predicts phase lifetimes for YSOs 
of ~ ICyrs, falling ofi' dramatically for objects with iboi^ 
lO^I/Q. The compact Hll-regions in the .RMS' survey have 
phase lifetimes between 2-4xl0^yrs. 

In the future, the framework of the simulations pre- 
sented here could be adapted to investigate Galactic struc- 
ture, once the distances to the RMS survey objects are 
known to a greater precision. Future modelling of the RMS 
population will take into account the evolution of the SED. 
This will include the increasing luminosity and heating as 
well as the dispersal of circumstellar material by feedback. 
Bipolar cavities in particular make the SEP dependent on 
the viewing angle jWhitnev et al.l |2003| ) and a simulation 
of a large sample with random viewing angles will be con- 
strained by comparison with the colour distribution as well 
as properties of the CO outflows and near-IR reflection neb- 
ulae. 
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Figure Al. The K-hand extinction of a sample of Galactic clus- 
ters versus the inferred gas column density from our model of the 
Galactic gas distribution. 
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APPENDIX A: COMPUTATION OF 
LINE-OF-SIGHT EXTINCTION 

In calculating the extinction towards each object in our sim- 
ulation, we cannot simply use 2-D extinction maps of the 
Galactic Plane, as these do not give any indication of how 
extinction varies with distance along a given sightline. In- 
stead, we use test points with known distances and extinc- 
tions, and derive a relation between the intervening gas col- 
umn density and the line of sight extinction by assuming 
that the extinction is some function of the interstellar gas 
column density between the Earth and the object. This rela- 
tion should be linear, as they are both directly proportional 
to the optical depth. 

For the test points, we must use objects that span a 
large range of distances and extinctions, as we will need to 
use this function to compute the extinctions to objects at 
the far side of the Galaxy. For this reason, we have chosen 
to use young stellar clusters as the test points. They are 
very bright in the near-IR, and so can be observed at very 
large extinctions (Ay ~ 30); while spectroscopic observa- 
tions of the stellar content allow the determination of both 
spectrophotometric and kinematic distances. 

In Table IXn we list a sample of clusters with known dis- 
tances and if-band extinctions which we use to derive the 
relation between extinction and column density. For each 
cluster, we use the model of the Galactic electron density 
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Figure A2. Illustration of the 50% completeness level of the MSX survey at 21/im, as a function of both I and b. 
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Table Al. Compilation of young Galactic clusters used in the cal- 
ibration of free-el ectron column dens ity with K-band ext i nction . 
R eferences are, 1: iBibby et al] ||200^: 2: iMessineo et al.l ||2008|): 
3:lMessineo et al.l ||2009^: 4:lDayies et al.l ll2007l'): 5: IPavies et al.l 
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described in Sect. 12.11 to calculate the integrated gas col- 
umn density along the line-of-sight from Earth to the clus- 
ter, . We also compleme nt these data with tha t from 
cepheids, using the data from iBerdnikov et al.l (|2000l ). 

In Fig. lAll we plot Ak against Eg- for each cluster in 
the sample. For the cepheids, which are all at low column 
densities, we take the average of the whole sample. As ex- 
pected, the relation is consistent with being linear, where 
we find, 



(-0.15 ± 0.08) + (1.24 ± 0.12) x (S^- /lO^ 



)(A1) 



Note that if the gas distribution is altered this relation 
must be recalculated. To convert the if-band extinction Ak 
to the extinction at 21/im A21, we must know the ratio 
A21/AK for interstellar extinction. Though not well stud- 
ied, there seems to be some agreement in the literature that 
this value is -^0.4-0.6 from studies of evolved stars, star- 



forming regions, t he Galactic Centre, and the diffuse in- 
terste llar medium | Drainc' 1989'; 'Lutz 1999'; Fl aherty et al.l 

.20071 ; IChapman et al 2009 ; .Messineo et al.. .20051 ). For the 

remainder of this work, we assume A21/AK = 0.5. 



APPENDIX B: COMPLETENESS OF THE MSX 
SURVEY AT 2lAtm 

In order to accurately predict whether or not a source would 
be picked up by the MSX survey, we must first know the 
detection limits of the survey. The survey is not isotropic, 
as some regions of the sky were scanned more times than 
others, so we must determine the completeness as a function 
of viewing angle. 

To measure the completeness of the MSX survey, 
we took representitive cutout images l°xl° in size. We 
then analyze d the images with the photometry package 
STARFINDER (|Diolaiti et al.l bOOd ). The input parameters 
were adjusted until the closest match was achieved between 
the STARFINDER and MSX source lists. The typical overlap 
between the sources was ^80%. 

Using a point-spread-function (PSF) measured from the 
test image, we then inserted fake stars into the image of a set 
brightness. No more than 40 stars were added at one time so 
as not to affect the crowding. We then ran the STARFINDER 
algorithm on the image and measured the recovery rate. The 
process was repeated five times per fake-star brightness in- 
terval, and six different brightness levels were tested between 
12Jy and 2Jy. This process was repeated on a large number 
of fields, spaced at intervals of 5° in / and 2° in 6, a total of 
345 images in all. 



